#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <common.h>
#include <ffdot.h>
#include <parse_config_file.h>
#include <parse_init_cond.h>
#include <evolve.h>
#include <plummer_model.h>

int main( int argc, char *argv[] ) {
  
  if( argc != 4 ) {
    fprintf( stderr, "Usage: ./nbodygpu <config file> <initial conditions> <output file>\n" );
    return 0;
  }
  
  SIM_DATA sim0;  
  
  FILE *in_config;
  in_config = fopen( argv[1], "r" );
  FILE *in_coords;
  in_coords = fopen( argv[2], "r" );
  FILE *out_coords;
  out_coords = fopen( argv[3], "w" );
  
  // Load the configuration file of the simulation
  parse_config_file( in_config, &sim0 );
  
  // Allocate neccesary memory
  allocate_arrays( &sim0 );
  
  // Check everything is in order
  check_simulation_params( &sim0 );
  
  // Load initial conditions
  parse_ic_file( in_coords, &sim0 );
  
  // Initializate timesteps using Aarseth's procedure (see nbody0.f)
  initialize_timesteps( &sim0 );
  
  // Evolution loop
  double tnext = 0.0;
  
  double E0 = energy_check( &sim0 ); 
  double E  = E0;
  double new_E = 0.;
  
  while( sim0.Time < sim0.tcrit ) {
    
    int L = build_blocksteps( &sim0 );
    
    predict( &sim0 );
    
    advance_block( &sim0, L );
    
    //Diagnosis
    if( sim0.Time > tnext ) {
      
      new_E = energy_check( &sim0 ); 
      fprintf( stderr, "TIME = %lf NSTEPS = %lf  dE0/E0 = %4.4e dE/E = %4.4e\n", 
                      sim0.Time, sim0.nsteps, fabs(new_E/E0-1.), fabs(new_E/E-1.) );
      
      E = new_E;
       
      dump_phasecoordinates( out_coords, &sim0 );    
      tnext += sim0.dtnext;
    }
  }
  
  fclose( in_config ); 
  fclose( in_coords ); 
  fclose( out_coords );
  
  return 0; 
}
